Introduction

A previous study had conducted both gene expression and metabolomics profiling of tissue samples from breast cancer patients. This vigenette will highlight the analysis we conduct on the breast cancer data.

Loading in IntLIM and Files

IntLIM, available from Github, can be installed as in the documentation. Once IntLIM is installed, what is necessary is loading in the library.

rm(list = ls())
library(IntLim)

For breast cancer study, both gene expression (available on Gene Expression Omnibus Accession Number: GSE27751) and metabolomics (http://www.jci.org/articles/view/71180/sd/2) data are available online. Much of this data had been processed as previously described (3). Probes from the Affymetrix data not mapping to a gene symbol were removed. Additionally, as with NCI-60 data, only the probe corresponding to the highest mean expression was used for analysis when multiple probes corresponded to a single gene. This resulted in a total of 20,254 genes for 108 patient samples. The Metabolon data did not need to be filtered by coefficient of variation, as there were no technical replicates. The resulting data consisted of 536 metabolites with 132 patient samples.

This data has been formatted for IntLIM. We load in the data as follows. The bc.ambs.csv meta file contains a list of phenotypic data file, metabolite data file, gene expression data file, metabolite meta file, and gene expression meta file. The function OutputStats will give a summary of the NCI-60 data

inputData <- IntLim::ReadData('bc.ambs.csv',metabid='id',geneid='id')
## Setting options('download.file.method.GEOquery'='auto')
## Setting options('GEOquery.inmemory.gpl'=FALSE)
## [1] "CreateMultiDataSet created"
IntLim::OutputStats(inputData)
##   Num_Genes Num_Metabolites Num_Samples_withExpression
## 1     20254             536                        108
##   Num_Samples_withExpression.1 Num_Samples_inCommon
## 1                          132                  108

From the OutputStats, we find that we have gene expression data involving 20,254 genes with 108 patient samples and metabolite abundance data involving 536 metabolites with 132 patient samples.

Filtering Gene Expression and Metabolite Data

We remove genes with mean belows below the 10th percentile. Furthermore, we remove metabolites with more than 80% missing values. This results in gene expression data involving 18,228 genes and 108 patient samples and metabolite abundance data involving 379 metabolites and 132 patient samples.

inputDatafilt <- IntLim::FilterData(inputData,geneperc=0.10, metabmiss = 0.80)
## [1] "No metabolite filtering by percentile is applied"
IntLim::OutputStats(inputDatafilt)
##   Num_Genes Num_Metabolites Num_Samples_withExpression
## 1     18228             379                        108
##   Num_Samples_withExpression.1 Num_Samples_inCommon
## 1                          132                  108

We can obtain boxplot distributions of the data as follows. This is used to make figures.

IntLim::PlotDistributions(inputDatafilt)

Principal Component Analysis

The principal component analysis is performed on filtered metabolite and gene expression data to obtain visual representations showing how different sub-sections of the data could be grouped into different clusters. Common samples patient samples (either tumor or adjacent non-tumor samples). Blue samples indicate tumor samples and red samples indicate non-tumor samples. Note the clear delineation of samples.

PlotPCA(inputDatafilt, stype = "DIAG", common = T)
## Warning in RColorBrewer::brewer.pal(numcateg, palette): minimal value for n is 3, returning requested palette with 3 different levels

Running Linear Model

The linear model is for integrating transcriptomic and metabolomics data is: E(m|g,t) = β1 + β2 g + β3 p + β4 (g:p) + ε where ‘m’ and ‘g’ are log-transformed metabolite abundances and gene levels respectively, ‘p’ is phenotype (cancer type, patient diagnosis, treatment group, etc), ‘(g:p)’ is the association between gene expression and phenotype, and ‘ε’ is the error term that is normally distributed. A statistically significant p value of the ‘(g:p)’ association term indicates that the slope relating gene expression and metabolite abundance is different from one phenotype compared to another. We run a linear model on tumor (n = 61) and non-tumor samples (n = 47) that included 18,228 genes and 379 metabolites (total of 6,908,412 possible associations and hence models). For genes and metabolites that had standard deviations of 0 in one of the groups, we assign a p value of NA. The model is run as below by calling RunIntLim(). DistPvalues() allows us to obtain a distribution of p-values for the (g:p) term.

myres <- IntLim::RunIntLim(inputDatafilt, stype="DIAG")
## [1] "Running the analysis on"
## 
## NORMAL  TUMOR 
##     47     61 
## [1] "10 % complete"
## [1] "20 % complete"
## [1] "30 % complete"
## [1] "40 % complete"
## [1] "50 % complete"
## [1] "60 % complete"
## [1] "70 % complete"
## [1] "80 % complete"
## [1] "90 % complete"
##    user  system elapsed 
## 347.855  33.232 382.035
IntLim::DistPvalues(myres)

The next step is to process the results of this model by filtering the results of the linear model by FDR-adjusted p-value cutoff (0.10 selected here) for the (g:p) association coefficient and calculate the correlations of the gene-metabolite pairs in each group (tumor and non-tumor) for the filtered results. We further may only interested in results that have an absolute correlation value difference (0.50 selected here). This is done with the ProcessResults function. In addition we also develop a heatmap of the gene-metabolite association correlations for the selected groups.

myres <- IntLim::ProcessResults(myres,  inputDatafilt, diffcorr = 0.50, pvalcutoff = 0.05)
IntLim::CorrHeatmap(myres)
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
dim(myres@corr)
## [1] 2842    4

From this model we find 2842 gene-metabolite correlations that have an association FDR-adjusted p-value of 0.05 and an absolute value correlation difference of 0.5 or greater. The top pairs are shown below.

corr.table <- myres@corr
abs.corrdiff <- abs(myres@corr$NORMAL - myres@corr$TUMOR)
sort.table <- corr.table[order(-abs.corrdiff),]
sort.table[1:10,]
##                                   metab     gene     NORMAL      TUMOR
## 445                   adrenate (22:4n6)   PDLIM4  0.5925755 -0.6262295
## 1035 docosapentaenoate (n3 DPA; 22:5n3) HIST1H3A -0.5074006  0.6666402
## 454                   adrenate (22:4n6) HIST1H3A -0.4935816  0.6429931
## 472                   adrenate (22:4n6)     GLI3  0.5703712 -0.5602327
## 334                   adrenate (22:4n6)    GRIA4  0.5279288 -0.6001586
## 1372                          guanosine   IGFBP4  0.5839669 -0.5426230
## 378                   adrenate (22:4n6)   IGFBP4  0.5797386 -0.5453728
## 1537                            leucine  COL14A1  0.5366559 -0.5864622
## 1032 docosapentaenoate (n3 DPA; 22:5n3)   PDLIM4  0.5124884 -0.6061952
## 416                   adrenate (22:4n6)     FHL2  0.5068810 -0.6087784

We can show some example plots of some of these pairs. The first example is the PDLIM4 vs. adrenate (22:4n6). There appears to be an error with the scatterplot.

IntLim::PlotGMPair(inputDatafilt, stype = "DIAG", geneName = "PDLIM4", metabName = "adrenate (22:4n6)")